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Abstract 

We develop an efficient numerical method to study the quantum critical behavior of disordered systems 
with 0(N) order-parameter symmetry in the large -A^ limit. It is based on the iterative solution of the 
large-A^ saddle-point equations combined with a fast algorithm for inverting the arising large sparse ran- 
dom matrices. As an example, we consider the superconductor-metal quantum phase transition in disor- 
dered nanowires. We study the behavior of various observables near the quantum phase transition. Our 
results agree with recent renormalization group predictions, i.e., the transition is governed by an infinite- 
randomness critical point, accompanied by quantum Griffiths singularities. Our method is highly efficient 
because the numerical effort for each iteration scales linearly with the system size. This allows us to study 
larger systems, with up to 1024 sites, than previous methods. We also discuss generalizations to higher 
dimensions and other systems including the itinerant antiferomagnetic transitions in disordered metals. 
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1. Introduction 

Randomness can have much more dramatic effects at quantum phase transitions than at classical phase 
transitions because quenched disorder is perfectly correlated in the imaginary time direction which needs 
to be included at quantum phase transitions. Imaginary time acts as an additional coordinate with infinite 
extension at absolute zero temperature. Therefore, the impurities and defects are effectively very large 
which leads to strong-disorder phenomena including power-law quantum Griffiths singularities 111 |2|, |3|], 
infinite-randomness critical points characterized by exponential scaling 13 Eh. ™d smeared phase transi- 
tions For example, the zero-temperature quantum phase transition in the random transverse-field Ising 
model is governed by an infinite-randomness critical point featuring slow activated (exponential) rather 
than power-law dynamical scaling. It is accompanied by quantum Griffiths singularities. This means, ob- 
servables are expected to be singular not just at criticality but in a whole parameter region near the critical 
point which is called the quantum Griffiths phase. 

Quantum Griffiths singularities are caused by rare spatial configurations of the disorder Due to statisti- 
cal fluctuations, one can always find spatial regions (rare regions) which are impurity free. The probability 
ViV^B.) to find such a rare region is exponentially small in its volume Vrr, !P(Vrr) ~ exp(-^yRR) with b 
being a constant that depends on the disorder strength. Close to a magnetic phase transition, the rare region 
can be locally in the magnetic phase while the bulk system is still non-magnetic. When the characteristic 
energy e of such a rare region decays exponentially with its volume, e ~ exp(-cyRR) (as in the case of the 
transverse-field Ising model), the resulting rare-region density of states has power-law form, p(e) ~ e''"', 
where A - b/c is the non-universal Griffiths exponent. A takes the value zero at the quantum critical point 
and increases throughout the quantum Griffiths phase. The singular density of states of the rare regions 
leads to quantum Griffiths singularities of several thermodynamic observables including order-parameter 
susceptibility, ~ T'* specific heat, C ~ T^, entropy, S ~ T'', and zero-temperature magnetization-field 
curve m ~ /? ' (for reviews see, e.g., Refs. JtIIMI)- 
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Many interesting models in statistical mechanics and field theory contain some integer-valued parame- 
ter and can be solved in the large-A^ limit. Therefore, the large-A^ method is a very useful tool to study 
classical and quantum phase transitions. An early example is the Berlin-Kac spherical model |@] which 
is equivalent to a classical 0(N) order parameter field theory in the large-A^ limit ifioll . Analogously, the 
quantum spherical model 111 Ull2lll3ll has been used to investigate quantum critical behavior In both cases, 
is the number of order parameter components. Another potential application of the large-A^ method are 
S^({N) Kondo models \IM with spin-degeneracy A^. In all of these cases, the partition function can be 
evaluated in saddle point approximation in the limit » 1, leading to self-consistent equations. In clean 
systems, these equations can often be solved analytically. However, in the presence of disorder, one obtains 
a large number of coupled self-consistent equations which can be solved only numerically. 

In this paper, we develop a new efficient numerical method to study critical behavior of disordered sys- 
tem with 0(N) order-parameter symmetry in the large -A^ limit. We apply this method to the superconductor- 
metal quantum phase transition in disordered nanowires. Using a strong-disorder renormalization group, it 
has recently been predicted that this transition is in the same universality class as the random transverse- 
field Ising model. We confirm these predictions numerically. We also find the behaviors of observables as 
a function of temperature and an external field. They follow the expected quantum Griffiths power laws. 
We consider up to 3000 disorder realizations for system sizes L = 256 and 1024. The paper is organized 
as follows: In Sec.|2]we introduce the model: a continuum Landau-Ginzburg-Wilson order-parameter field 
theory in the presence of dissipation; and we generalize the theory to quenched disordered systems. Then, 
we discuss the predicted critical behavior of this model and derive the large-A^ formulation. In Sec. |3] we 
review an existing numerical approach to this model. In Sec. [H we present our numerical method to study 
the quantum critical behavior. We discuss the results in Sec.|5] and we compare them to the behavior pre- 
dicted by the strong-disorder renormalization group. Sec. |6]is devoted to the computational performance 
of our method. Finally, we conclude in Sec. |2]by discussing and comparing our numerical method to the 
existing one. We also discuss generalizations to higher dimensions and other models. 



2. The model 

We start from the quantum Landau-Ginzburg-Wilson free-energy functional for an A'^-component vec- 
tor order parameter (p in one space dimension. For a clean system with overdamped order parameter 
dynamics the Landau-Ginzburg-Wilson action reads|3 



I r r ' u 

yT ^ r 2 r r'^^ 

+ ^ l<j^nl J t/x|^(x, w„)| —h J dx J dT(p(x,T), 



(1) 



where a is the bare distance from criticality. y and J are the strength of dissipation and interaction, re- 
spectively, u is the standard quartic coefficient, h is a uniform external field conjugate to the order param- 
eter ifi{x, (x>„) is the Fourier transform of the order parameter (p{x, r) with respect to imaginary time, and 
u)n - 2miT is a Matsubara frequency. The above action with N - 2 order parameter components (equiv- 
alent to one complex order parameter) has been used to describe 11511 the superconductor-metal transition 
in nanowires lll6ll . This transition is driven by pair-braking interactions, possibly due to random magnetic 
moments trapped on the wire surface lll6ll . which also introduce quenched disorder in the nanowire. The 
action ([T]i can be generalized to d = 3 space dimensions and N - 3 order parameter components, in this 
case, it describes itinerant antiferromagnetic quantum phase transitions iITtIIisII . 

In the presence of quenched disorder, the functional form of Eq. ^ does not change qualitatively. 
However, the coupling constants become random functions of position x. The full effect of disorder can 
be realized by setting m = y = 1 while considering the couplings a and J to be randomly distributed in 
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space ifigl] . The quantum phase transition in zero external field can be tuned by changing the mean of the 
a, distribution, a. 

Recently, the model ([1]) has been investigated by means of a strong-disorder renormahzation group 
method i20l Bill . This theory predicts that the model falls in the same universality class as the one- 
dimensional random transverse-field Ising model which was studied extensively by Fisher Jstl- Thus, the 
phase transition is characterized by an infinite-randomness critical point at which the dynamical scaling is 
exponential instead of power-law. Off criticality, the behaviors of observables are characterized by strong 
quantum Griffiths singularities. 

Let us focus on the Griffiths phase on the disordered side of the transition, where the distance from 
criticality 6 = a - cic > 0- The strong-disorder renormahzation group predicts the disorder averaged 
equal-time correlation function C{x) to behave as ^ 

^, . exp[-(x/g)-(27;rV4)'/3(x/^)'/3] 

C(x) — (2) 

for large distances x. Here, ^ is the correlation length which diverges as ^ ~ \6Y^ with v = 2 as the critical 
point is approached. The disorder averaged order parameter as a function of the external field h in the 
Griffiths phase has the singular form [53 

^{h) ~ . (3) 

Here, A is the non-universal Griffiths exponent which vanishes at criticality as A ~ 6^''^ with critical exponent 
i// - 1/2. Right at criticality, the theory predicts logarithmic behavior rather than a power law fl, 

<p(h)~[log{ho/h)r-'''^. (4) 

Here, the exponent 4> — (I + V5)/2 equals to the golden mean, and ho is some microscopic energy scale. 

The average order parameter susceptibility as a function of temperature T in the disordered Griffiths 
phase is expected to have the form 

x(T) ~ r-"-' (5) 

with the same /l-exponent as in Eq. (|3]l. 

Our goal is to test the strong-disorder renormahzation group predictions by means of a numerical 
method. As a first step, we discretize the continuum model ([TJ in space and Fourier-transform from imagi- 
nary time T to Matsubara frequency lo„. The discretized Landau-Ginzburg-Wilson action has the form 

/■= 1 to,, 
i=l ojii 

where L is the system size. The nearest-neighbor interactions /, > and the mass terms ff, (bare local 
distances from criticahty) are random quantities. The critical behavior of the model (|6]l can be studied 
in the limit of a large number of order parameter components A^. In this limit, the above action can be 
reduced to a Gaussian form. This can be done in several ways, for example by decomposing the square 
of each component of the order parameter \ipf\tL)„)f- into its average {(p^) and fluctuation A\ipf\a>„)f': 
|<jo':*'(w„)p = {(p~} + Al^'^^w,,)!^. Substituting this into the quartic term of the action (|6) and using the 
central limit theorem, the quartic term can be replaced by {i^-)|ijo,(w„)p. This leads to the Gaussian action 
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i,j=l ui„ i=\ 



The coupling matrix is given by 



Mij = -JiSij+x - Jj6ij-i + (n + 2Ji)6ij . (8) 
The renormalized local distance r, from criticality at site / must be determined self-consistently from 



n = ai + {ifi") , (9) 



where {ip^) is given by 



L 

{ip") = r ^[M + |^^„|i],i + Yj ^u^ik ■ (10) 

w„ j,k=\ 

Here, 1 is the identity matrix. In the presence of disorder, the self -consistent equations (|9]i at different 
sites are not identical. We thus arrive at a large number of coupled non-linear self-consistent equations. 
Therefore, numerical techniques are required to solve them. 

3. Existing numerical approach 

In this section, we review the numerical method proposed by Del Maestro et al. to study the 
model (|7]i at zero temperature and in the absence of an external field {h = 0). The matrix M is spectral 
decomposed in terms of its orthogonal eigenvectors Vij and eigenvalues e, as 

L 

Y^MijVjk^Vikek. (11) 
Using this decomposition, the inverse matrix in Eq. (fTOl i can be written as 

[M+Kii],' = y^^. (12) 



At zero temperature the sum over Matsubara frequencies in Eq. ( fTOt turns into an integral which can 
be performed analytically. This leads to the self-consistent equations (for h - 0), 

i jjyijf log (i + ^) + = • (13) 

Here, for convergence of the frequency integral, an ultra violet cutoff" A,^ is introduced. Numerical solutions 
to Eq. (fT3] ) were obtained by an iteration process using a modified Powell's hybrid method. The method 
works well for large distances from criticality and small system sizes, but it becomes computationally 
prohibitive near criticality where the correlation length ^ becomes of order of the system size. This problem 
can be partially overcome by implementing a clever iterative solve-join-patch procedure. However, the 
system size L is still limited because large matrices need to be fully diagonalized which requires 0(L?) 
operations per iteration. Therefore, for large L the method gets very slow. 



As the result, the largest sizes studied in Ref. 112211 were L = 128. The authors analyzed equal time 



correlations, energy gap statistics and dynamical susceptibilities and found them in agreement with the 



strong-disorder renormalization group predictions [20 21]. The method was also used in Ref. |23l to study 
the conductivity. 
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4. Method 



We now present a novel numerical method to study the model ([7) at non-zero temperatures. Its numer- 
ical effort scales linearly with system size L (per iteration) compared with the scaling of the numerical 
method outlined in Sec. [3] The basic idea of our method is that, for h - 0, we only need the diagonal 
elements of the inverse matrix M ' to iterate the self-consistent Eq. (|9]l- The numerical effort for finding 
the diagonal elements of the inverse of a sparse matrix is much smaller than that of a full diagonalization. 
Combining Eqs. (|9]l and ( fTOt . the system of self-consistent equations at non-zero temperatures T, and in 
the presence of an external field h, reads 



r,- = 2T YjiM + 2TmT\\f + TU'^ + /z^ ^ + • (14) 

Here, m - h^{2jiTY'^ with an ultra-violet cutoff frequency A„. To solve these equations ( fT4l i iteratively, 
we find the inverses of the tridiagonajl matrices [M,; + InnTl] and M,j using the fast method proposed in 



Ref. 112411 . This algorithm is summarized in Appendix A In zero external field, we only need the diagonal 



elements of [M^ -i- ImiTty and the number of operations per iteration scales linearly with system size 
L, while it scales quadratically in the presence of a field because for h + 0, full inversion of the matrix is 
required. 

Once the full set of r-, has been obtained, we can compute observables from the quadratic action Q. 
Let us first consider observables in the absence of an external field. The equal-time correlation function 
C{x) - {(f x{T)(p[{T)) averaged over disorder realizations can be obtained from Eq . (|7]i. 



C(x) = 



L-X { m 

-V 



^2[M-i-27r«ri]-iy+, + 



M7 



(15) 



where the overbar indicates the average over disorder configurations. Similarly, in the zero external field, 
we can calculate the order parameter susceptibility as a function of temperature. The disorder-averaged 
order parameter susceptibility x^T^) can be expressed as 



(16) 



1=1 k=\ 



In the presence of an external field, we need to include h in the solution of Eq. (fT4l i. We can then 
compute the order-parameter vs. field curve. The disorder-averaged order parameter reads 



(17) 



1=1 k=l 



We note that the number of operations to calculate observables for one disorder configuration scales 
quadratically with the system size L. However, this needs to be done only once, outside the loop that iterates 
the self-consistent equations. At low temperatures, according to Eq. (fT4l i. we need to invert a huge number 
of matrices [M,j + ImiTl] per iteration (one for each Matsubara frequency). Naively, one might therefore 
expect the numerical effort to scale linearly in 1/T. However, these matrices are not very different. We can 



therefore accelerate the method by combining them appropriately. This is explained in Appendix B 



We use open boundary conditions. 
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Figure 1 : (Color online) Order-parameter susceptibility x versus temperature T for various distances from ciiticality 5- in the Griffiths 
phase. All data are averaged over 3000 disorder configurations with system size L = 256. The solid lines represent fits to the Griffiths 
power law l|5), x(X) ~ ' • °ver the temperature range T = 10"^ - 1 .5 X 10"^. 



5. Results 

In this section, we report results of our numerical calculations of the model O. We consider the 
interactions 7, to be uniformly distributed on (0, 1) with mean J - 0.5 and the bare local distances from 
criticality a, to be Gaussian distributed with mean a and variance 0.25. 

An advantage of our method is that it gives direct access to the temperature dependencies of observ- 
ables. For example, we calculate the zero-field order parameter susceptibility as a function of temperature 
for various values of the control parameter a according to Eq. ( fT6l l. At low temperatures, the Griffiths 
power law (|5]l describes the data very well (see Figure [!}. The non-universal Griffiths exponent A can be 
determined from fits in the temperature range T - 10^^ - 1 .5 x 10"~. Figure|2a) shows how A varies as the 
distance from criticality 6 = a-ac changes. The power law A ~ S^"^ describes the data well with the critical 
point Uc - -0.85(3), and exponents v = 2.0(2) and i// = 0.51(2). Here, the number in brackets indicates the 
uncertainty in the last digit. These results are consistent with the predictions of Refs. |2^ 21 ]. 



We also compute the order parameter as a function of an external field T — 10 for various a 
(Figure|3]l. The off-critical data {6 > 0) are described by the Griffiths power law (|3]l with an exponent A. At 
the critical point, the (p{li) curve follows the logarithmic dependence (|4|i with exponents iff - 0.51(2) and 
- 1.61(2). The value for exponent (p is in agreement with the predicted one lEolElll . The values of the 
Griffiths exponent A match those extracted from susceptibility data (see Figure |2] (a)). The deviation near 
the critical point may be due to the fact that the correlation length becomes comparable to the system size 
and correspondingly causes finite-size effects in the data. 

In addition, in the absence of an external field h, for system size L - 1024, we compute the disorder- 
averaged correlation functions (flSt at temperature T = 10 for various values of a (see Figure |4]i. The 
values of correlation length ^ can be extracted by fitting the data to Eq. We find good agreement of 
the data with Eq. ^ for distances between x = 5 and some cutoff at which the curves start to deviate from 
the zero-temperature behaviors due to temperature effects and where curves start to become noisy because 
correlations become dominated by very rare large clusters. 

Figure|2tb) shows how the correlation length ^ changes with distance from criticality 5. The data can be 
fitted to the power law ^ ~ |(5| as expected 151]. By fitting, we extract the critical point Uc - -0.85(3) and 
exponent v = 2.0(2). The values of exponent v and critical point are in agreement with those obtained 
from;^f(r) and tp(h). 
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Figure 2: (Color online) a) The Griffiths exponent A versus distance from criticality S. The solid line is a fit to the power law A ~ S'^^'. 
b) The correlation length f obtained by analyzing correlation function data versus distance 6 from criticality. The solid line is a fit to 
a power law, resulting in a ciitical point of = -0.85(3) and the correlation length exponent v = 2.0(2). 



6. Computational performance 

In this section, we discuss the execution time of our method for solving the self-consistent Eqs. (fT4l l 
iteratively (i.e., the time needed to get a full set of renormalized distances from criticality r,). In our method, 
the time per iteration scales linearly with the system size L in the absence of an external field because 
the operation count is dominated by the matrix inversion. Thus, the disorder-averaged execution time 
f ~ «itL for a single disorder configuration, where njt is the number of iterations needed for convergence 
of the self-consistent Eqs. ( fT4l i. The number of iterations nit depends on the disorder configuration, it is 
larger for a disorder realization which has locally ordered rare regions with smaller a. In the conventional 
paramagnetic phase, i.e., for larger values of a away from criticality, locally ordered rare regions are almost 
absent, therefore the number of iterations «it is a constant. Thus, in the conventional paramagnetic phase, 
the execution time is expected to scale linearly with the system size, t ~ L. Figure |5] shows that it indeed 
scales linearly with the system size for a = 1 . In contrast, in the quantum Griffiths phase, where locally 
ordered rare regions are present, rin is expected to be large and to become larger close to criticality. If we 
compare two different system sizes in the quantum Griffiths phase, the larger system is expected to have 
locally ordered rare region with higher probability. Thus, in the quantum Griffiths phase the number of 
iterations «it is expected to be a function of system size L, which we model as «it ~ U with some non- 
negative exponent y. Therefore, in the quantum Griffiths phase the execution time does not scale linearly 
with the system size but it behaves as r ~ L''^'. Figure |5] shows that for a - -0.6 in the quantum Griffihts 
phase, the disorder averaged execution time t does not scales linearly with L but behaves as power law 
f ~ L'+i withy = 0.6. 

Because our method performs the Matsubara sums numerically, the effort increases with decreasing 
temperature T. As shown in [Appendix B[ this increase is only logarithmic in l/T if we approximately 
combine higher Matsubara frequencies. 



7. Conclusions 

In summary, we have developed an efficient numerical method for studying quantum phase transitions 
in disordered systems with 0{N) order parameter symmetry in the large -A^ limit. Our algorithm solves it- 
eratively the large- self-consistent equations for the renormalized distances from criticality using the fast 
method of Ref. 12411 for the necessary matrix inversions. We have applied our method to the superconductor- 
metal quantum phase transition in nanowires and studied the critical behavior of various observables near 



7 



-0.85 

■ Fit to Eq.(3) 

■ Fit to Eq.(4) 



Figure 3: (Color online) Order parameter if) versus external field h for various a. The data ai'e averaged over 3000 disorder config- 
urations of system size L = 256. In the field range h = 10"'' to 2 X 10"', the dotted and solid lines represent fits to Eq. l|4) and the 
Griffiths power law l|3], respectively. 



the transition. Our results are in agreement with strong-disorder renormalization predictions H El that 
the quantum phase transition is governed by infinite-randomness critical point accompanied by quantum 
Griffiths singularities. 

Let us compare the performance of our method with that of the method proposed in Ref. ll22ll and out- 
lined in Sec. |3] The main difference is how the sums over the Matsubara frequencies in the self-consistent 
equations (|9j are handled. The method of Ref. ll22ll works at T = where the Matsubara sum becomes 
an integral. This integral is performed analytically which saves computation time. However, the price is a 
complete diagonalization of the coupling matrix M which is very costly {0{Lr') operations per iteration). 
Moreover, observables at T it are not directly accessible. 

In contrast, our method performs the Matsubara sum numerically which allows us to use the fast matrix 
inversion of Ref. ll24ll (which needs just 0{L) operations per iteration) instead of a full diagonalization. 
Furthermore, we can calculate observables at T ^t 0. However, our effort increases with decreasing T . 
Thus, the two methods are in some sense complimentary. The method of Ref. [22] is favourable for small 
systems when true T - Q results are desired. Our method works better for larger systems at moderately 
low temperatures. 

We also emphasize that all our results have been obtained by converging the self-consistent equations 
(|9]l by means of a simple mixing scheme. Even better performance could be obtained by combining our 
matrix inversion scheme with the solve-join-patch algorithm ll22ll for convergence acceleration. 

Our method can be generalized to higher-dimensional problems. The self-consistent equations can be 
solved in the same way, using a fast method for inverting the arising sparse matrices. For two dimensional 



is a total number of sites. We therefore expect the cost of our method to scale as N^'^^^^ or A^^^^ in the 
quantum Griffiths and quantum paramagnetic phases, respectively. For three dimensional systems, sparse 
matrices can be inverted in O(N^) operations \ 2^, correspondingly the cost of our method is expected to 
behave as A^''^^ (A^ is number of sites) in the quantum Griffiths phase. In the quantum paramagnetic phase 
it should scale as A^-. 

A possible application of our method in three dimensions is the disordered itinerant antiferromagnetic 
quantum phase transitions lEolElll . The clean transition is described by a Landau-Ginzburg- Wilson theory 
which is generalization of the action ^to d - 3 space dimensions and N = 3 order parameter components 



|17|,[18|]. Introducing disorder leads to random mass terms as in the case of the superconductor-metal 



quantum phase transition in nanowires. 
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Figure 4: (Color online) The equal-time coiTelation functions for several values of a. All data are averaged over 3000 samples of size 
L = 1024 at T = 10"^ The solid lines are fits to Eq. (2)- Inset: Deviations of congelation function at fixed value of cc = -0.7 due 
to temperature effects and statistical error of an average over disorder configurations. The data represented by circles and stars are 
averaged over the same 1000 disorder configurations at T = 0.0025 and T = 10^^, respectively. The curves represented by triangles 
are averaged over different set of 1000 disorder configurations at 7 = 10"'. 



8. Acknowledgements 

This work has been supported by the NSF under Grant Nos. DMR-0906566 and DMR-1205803. 



Appendix A. Inversion of tridiagonal matrix 



In this Appendix we sketch the fast method for the inversion of a tridiagonal matrix outhned in Ref . Il24ll . 
The cost of finding the diagonal elements of the inverse matrix is 0{L) operations while inverting the full 
matrix costs 0{L^) operations. The basic idea is that the inverse matrix of the tridiagonal matrix Mjj can be 
represented by two sets of vectors vj and uj: Mj^ - UiVj. Let diagonal and ofFdiagonal elements of matrix 
Mij be Ma - a, and M, ,+i = M,;,+i = respectively. By combining a UL decomposition of the linear 
system for v and a UL decomposition of Mjj, one can determine the set of vectors 

1 b2---bi 

vi = — , ^i^-} 1 — T' ! = 2, (A.l) 

di di ■ ■ ■ di-\di 

where 

dn — a„ , di — flj — -p— , / = « - 1 , ■ ■ ■ , 1 . (A. 2) 

The set of vectors t/y can be found by combining a UL decomposition of the linear system for u and a 
UL decomposition of Mij, yielding 

1 b„-i+i ■■ - bn . /A ON 

Un = , Un-i = , !=1, (A.3) 

OnVn On-i ■ ■ ■ OnV,, 

where 

^/ 

6\ = a\ , 6i = ai '— , i — 2, - ■ ■ ,n . (A.4) 
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Figure 5: (Color online) At the temperature T = 10"' and in the zero field h = 0, execution time for a single disorder configuration t 
versus system size L for a = -0.6 and (5 = 1. All data are averaged over 1000 disorder realizations. The solid lines represent fits to 
the power-law. (times measured on an Intel Core 15 CPU) 

Finding both sets of vectors needs 0{L) operations, consequently the number of operations to extract 
the diagonal elements M^^' = m, v, of inverse matrix scales linearly with L while the cost of finding the full 
inverse matrix Mr.' = m/V; is 0{Lr). 

Appendix B. Acceleration of method 

In this Appendix we propose an approach to accelerate the summation over the Matsubara frequencies 
in our method. The idea is based on the fact that the critical behaviors are dominated by low-frequencies, 
correspondingly only matrices associated with low Matsubara frequencies <y„ have dominant contributions 
in Eq. (fT4l i. At higher w„, consecutive matrices change very little. Therefore, instead of finding diago- 
nal elements of [M^ + InTtilY^ for each Matsubara frequencies w,,, we invert matrices corresponding to 
n = 1, 100 and correspondingly calculating the sum of first 100 terms in Eq. (fT4l i exactly. Then, we 
approximate sum of the remaining terms corresponding to n > 100 (higher Matsubara frequencies) in the 
following way: we find diagonal elements of [M,-y + InTnl]"^ corresponding to the midpoints of subinter- 
vales obtained by dividing interval n - lO'^' -i- 1, lO'^^ (/ = 1, logio(OT/100)) into 90 subintervales 
of width lO'. Then, we approximate appropriate sum in Eq. (Il4l by summing over terms calculated at 
midpoints multiplied by lO'. In this case, numerical effort scales logarithmically as log[o(l/r) compared 
with l/T scaling in the case of exact summation. To check the magnitude of errors arising due to this 
approximation, we have compared observables calculated exactly and using acceleration method for the 
system with size L - 256 and control parameter a,- = -0.6 at the temperature T - 10"^. We have found 
that arising errors are less than 0.1%. 
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